Studying Soft Matter with "Soft" Potentials: Fast Lattice Monte Carlo Simulations 
and Corresponding Lattice Self-Consistent Field Calculations 
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The basic idea of fast Monte Carlo (MC) simulations is to perform particle-based MC simulations 
with the excluded-volume interactions modeled by "soft" repulsive potentials that allow particle 
overlapping. This gives much faster system relaxation and better sampling of the configurational 
space than conventional molecular simulations with "hard" repulsions that prevent particle over- 
lapping. Here we present fast lattice MC (FLMC) simulations for confined homopolymers, where 
multiple occupancy of lattice sites is allowed with a proper Boltzmann weight and thus the evalua- 
tion of nearest-neighbor interactions can also be avoided. When compared with the corresponding 
lattice field theories based on the same Hamiltonian, FLMC simulations further provide a powerful 
means for unambiguously and quantitatively revealing the correlation/fluctuation effects. 

PACS numbers: 05.10.Ln, 61.25.H-, 64.60.De 



Significant advances have been made in the past 
three decades in applying conventional, "particle-based" 
molecular simulation techniques (e.g., Monte Carlo (MC) 
and molecular dynamics simulations) to the study of 
polymeric systems. Even with coarse-grained models 
where each segment represents several to tens of re- 
peat units of real polymers, however, conventional molec- 
ular simulations of multi-chain systems are still hin- 
dered by "hard" excluded-volume interactions and, for 
off-lattice simulations, expensive pair-potential calcula- 
tions. The former is implemented in off-lattice simula- 
tions commonly by either hard-sphere repulsion or the 
harsh Lennard-Jones repulsion, and in lattice simula- 
tions by the self- and mutual-avoiding walk (SMAW). 
While such hard-core repulsions are needed to produce 
realistic dynamics, they greatly reduce the chain relax- 
ation towards equilibrium configurations as well as the 
efficiency of sampling the configurational space. On the 
other hand, the pair-potential calculations are expensive 
for dense polymeric systems, and lattice simulations are 
in general much faster and thus widely used. 

The basic idea of fast MC simulations is to perform 
particle-based MC simulations with the excluded-volume 
interactions modeled by "soft" repulsive potentials that 
allow particle overlapping (i.e., by soft particles whose in- 
teraction energy u(r) is finite when they overlap). This 
avoids the hard repulsions (i.e., u(r — ► 0) — ► oo as in the 
Lennard-Jones potential) used in conventional molecu- 
lar simulations, thus allowing much faster chain relax- 
ation and better sampling of the configurational space. 
Such coarse-grained models are particularly suitable for 
the study of equilibrium properties of soft matter such 
as polymers. It turns out that soft potentials are com- 
monly used in polymer field theories (e.g., the widely 
applied self-consistent field (SCF) theory with great 



'Electronic address: q.wang@colostate.edu 



success for many polymeric systems) [lj , where individ- 
ual polymer segments are modeled as volumeless points 
with the excluded-volume interactions described by ei- 
ther the Helfand compressibility L 2] or the incompress- 
ibility constraint. Using the same Hamiltonian in both 
fast MC simulations and polymer field theories therefore 
enables quantitative comparisons between them without 
any parameter-fitting to unambiguously reveal the con- 
sequences of approximations in the theories (e.g., the ef- 
fects of correlations and fluctuations neglected in SCF 
theory) 3. Soft potentials have also been used in dissi- 
pative particle dynamics [H and the dynamic mean-field 
density functional method || that focus on the system 
dynamics. 

In our recent wor k@|, fast off-lattice MC (FOMC) sim- 
ulations were performed for several systems, using an 
isotropic and position-independent soft repulsive poten- 
tial. As studied in detail there, while a spatial dis- 
cretization scheme used in previous FOMC simulations 
avoids pair-potential calculations, it is equivalent to an 
anisotropic and position- dependent pair potential that 
cannot be implemented in field theories. In this Let- 
ter, we report fast lattice MC (FLMC) simulations 
where multiple occupancy of lattice sites is allowed with 
a proper Boltzmann weight. While this Domb- Joyce 
modelQ was already studied in the early single-chain 
simulations of variable excluded-volume effects [9j], we 
demonstrate here its great advantages for multi-chain 
systems and further compare our FLMC results with the 
corresponding lattice self-consistent field (LSCF) calcu- 
lations based on the same Hamiltonian. 



As a simple example, we consider an inhomogeneous 
system of n homopolymer chains each consisting of N = 
60 segments confined between two parallel and impene- 
trable walls placed in the x direction. The energetic wall- 
polymer interaction is ignored, and the system Hamilto- 
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nian due to non-bonded interactions is given by 

(1) 
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where (3 = 1/ksT with ks being the Boltzmann constant 
and T the absolute temperature, the average segmental 
density po = nN/V with V being the total number of 
lattice sites, the microscopic density of polymer segments 
at a lattice site r is p(r) = Y^Jk=i S s =i ^r,R fe , s with Hk, s 
denoting the lattice position of the s th segment on the 
fc th chain and <5 r ,R fc „ being the Kronecker <5-function, and 
k = 0.48 denotes the compressibility. The value of N/k 
indicates that our system is nearly incompressible. Note 
that Eq. ([1]) can equivalently represent homopolymers in 
an implicit, good solvent [f| [HI • Hereafter we define two 
physical parameters that determine the system behavior: 
chain number density C = n/V and excluded- volume 
parameter B = N/kC. 

Our FLMC simulations are performed in a canoni- 
cal ensemble on a simple cubic lattice of V = L X L 2 
sites, where L x = 10 is the number of lattice sites in 
the x direction that can be occupied by polymers, and 
L = 20 ~ 60 (depending on the finite-size effects) is the 
number of lattice sites in both the y and z directions 
along which the periodic boundary conditions are ap- 
plied. Our trial moves include local moves (end-rotation 
and kink-jump), reptation, and bond-rotation (where a 
randomly selected bond vector of a chain is altered), oc- 
curring with probabilities of 0.1, 0.45, and 0.45, respec- 
tively. About 2 x 10 8 ~ 1.3 x 10 10 trial moves are per- 
formed in each simulation (depending on n and V), and 
the Metropolis acceptance criterion is used. We have 
estimated the error bar of each ensemble-averaged quan- 
tity in FLMC simulations as three times its standard de- 
viation, with the statistical correlations among samples 
collected after equilibration taken into account through 
the block analysis [ll|; the error bar is smaller than the 
symbol size in all the cases and thus not shown. 

To formulate the corresponding LSCF theory, we start 
from the canonical-ensemble partition function 

z = -. V eM-(3n c [{n k }]-pn E } (2) 

{Rfc, s } 

where the summation is over all possible lattice posi- 
tions of all polymer segments, and 7i is the Hamilto- 
nian due to chain connectivity; (3Ti c = if the chain 
connectivity is maintained for all chains on the lat- 
tice, and oo otherwise. After the Hubbard-Stratonovich 
transformation, i.e., inserting in Eq. ([2]) the identity 
1 = J ^<f>S>w exp{^ r u;(r)[po</! > (r) - /3(r)]} where <fi(r) 
is the (normalized) segmental density field constrained 
to p(r)/po and u)(r) is the (purely imaginary) conju- 
gate field imposing the constraint, we finally have Z — 
J @<f>!0Lie-xp{-nPf c [<p,uj]} with 
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E[0(r)-l] 2 -iEc(r)0(r)-lng[H (3) 



and the single-chain partition func- 
tion Q = J2{tl s } ex P[ — (3H c /n — 
£f =1 w(R s )/A r ]/£ { R s} exp(-/3H c /n), where we 

have omitted a constant factor in Z and re-scaled vari- 
ables according to iVcj(r) — » uj(r). The SCF equations 
are obtained by setting 6Pf c /5(j)(r) = 5f3f c /Suj(r) = 
(i.e., the mean-field approximation) and given by 



w(r) = (N/ K )[<t>(r) - 1} 
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exp [uj{r)/N] ^ 
( r ) = ^771 > , <ls(r)qN+i~s(r) 
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(4) 
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where q s (r) is the probability of finding a partial chain of 
s segments starting anywhere in the system and ending 
at r. According to the chain connectivity, we have 

q s +1 (r) = eXP[ - Uj{r)/N] ^s(r n ), q0 (r) = l (6) 



where r„ denotes a nearest neighbor of r on the lattice 
and z = 6 is the lattice coordination number. Finally, 



Q 



W (r) 



N 



q s (r)q N+1 - s {r) 



(7) 



Note that our LSCF formalism is different from that 
of Scheutjens and Fleer [l2|: In addition to the incom- 
pressibility constraint, their model uses nearest-neighbor 
interactions between different species. Since multiple oc- 
cupancy of lattice sites is allowed in our model, however, 
we can use S r . r > instead (refer to Eq. (Q])). This further 
avoids the evaluation of nearest-neighbor interactions in 
our FLMC simulation, which makes it very fast. 

For the confined homopolymers, we solve the SCF 
equations in ID using the Broyden method combined 
with a globally convergent strategy [l3j], where the resid- 
ual errors of Eq. ((4j at all x are less than 10 -12 . We 
then calculate the propagator Q PtS (x\x'), which corre- 
sponds to the probability of finding a partial chain of 
s + 1 segments that starts at x 1 and ends at x in the ob- 
tained conjugate field uj(x), from the following equation 
analogous to Eq. ([6]) 

Q PtS+ i(x\x') = cxp[-uj(x)/N]{z a Q PtS (x\x') 

+zi[Q PjS (x + l\x') + Q PtS (x - l\x')]} 

with the initial condition of Q p ^ s =o{x\x') — 8 X . X / and the 
appropriate boundary conditions, where zq = 2/3 is the 
fraction of nearest-neighbor lattice sites at the same x 
and Z\ = (1 — zq)/2. The mean-square chain radius of 
gyration in the x direction can finally be computed as 
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i r ,i '- x^(n „\ ^xJ2 X ' Qp,s( x \ x ')( x - x ') 2 



where the last term in the summation (i.e., with s = 
N — 1) is the mean-square chain end-to-end distance in 
the x direction. 
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FIG. 1: (Color online) Polymer segmental density profiles 
4>{x) from FLMC simulations at different chain number den- 
sities C, and from LSCF calculations (which correspond to 
C — > oo). x — 1 is the closest lattice layer to a confining wall, 
and x = 5 is in the middle of the film. 
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FIG. 2: (Color online) Ensemble-averaged non-bonded inter- 
action energy per chain (/3H E )/n and number of overlapping 
segment pairs occupying the same lattice sites n pa i r . 

Fig. [1] shows the ensemble-averaged segmental density 
profiles 4>(x) = (p(x))/po from FLMC simulations at dif- 
ferent chain number densities C; since 4>{x) is symmetric 
about the mid-plane between the two confining walls, it is 
shown only from x = 1 (closest to a wall) to 5 (in the mid- 
dle of the film). The smallest C value in FLMC simula- 
tions is obtained from a single chain with L = N, which is 
effectively self-avoiding due to the large corresponding B 
value. As C increases, 4>{x) increases near the walls and 
decreases in the interior of the film. This is in goo d agree- 
ment with recent field-theoretic simulations [10|. <f>(x) 
from LSCF calculations is also shown in Fig. [T] which 
corresponds to C — > oo and is approached by FLMC 
simulations with increasing C. For C > 0.1, 4>(x) is oscil- 
latory (i.e, exhibits local maxima at x = 2 and 4) rather 
than monotonic with x due to the packing effects. 

Fig. [5] shows how the ensemble-averaged non-bonded 



FIG. 3: (Color online) Log-log plot of mean-square chain radii 
of gyration (Rg,i) along different directions i as a function of 
chain number density C. 

interaction energy per chain ((3H E )/n varies with C. Its 
behavior can be well understood by considering the num- 
ber of overlapping segment pairs occupying the same 
lattice sites n pa ir also shown in the figure; note that 
f3H E /n = (N/2K)[(2/N 2 )(n paiy: /nC) + 1/NC - 1]. At 
small C values (< 0.002) we effectively have SMAW (i.e., 
(^pair) = 0) due to the large corresponding B values, 
while at large C values (n pa i r ) jnC approaches that pre- 
dicted by LSCF theory, which is obtained by equating 
j3H jn to the first term on the right-hand-side of Eq. ([3]) 
after LSCF equations are solved. Due to its mean-field 
approximation, LSCF theory underestimates n pa i r /nC, 
and gives n pair < for C < 1.664 x 10~ 2 . 

Fig. [3] shows how the mean-square chain radii of gy- 
ration (Rg^) (i = x,y, z) change from the (self-avoiding) 
single-chain case to the LSCF limit as C increases. In 
the directions parallel to the confining walls (i.e., i = y 
and z), (Rg i) monotonically decreases with increasing 
C and approaches the LSCF limit corresponding to the 
random walk where (R 2 gA ) = (N 2 - 1)/18N. In the direc- 
tion perpendicular to the walls, however, (R 2 x ) first in- 
creases with increasing C to a maximum located around 
C = 0.004, then decreases to approach the LSCF limit. 
The same behavior is also found in the mean-square chain 
end-to-end distances (data not shown). 

To summarize, we have presented for the first time 
FLMC simulation data ranging from the single-chain case 
all the way to the LSCF limit. In terms of the invariant 
degree of polymerization Af = [n/(V/Rg )] 2 , where the 
root-mean-square end-to-end distance of ideal chains is 
given by i? e , = VN-1 for our lattice polymers, these 
data correspond to M w 1.585 x 10~ 4 - 5.134 x 10 6 , 
more than ten orders of magnitude across. Even bet- 
ter, the total simulation time for obtaining these data is 
less than a couple of days on an Intel Core 2 Q6600 2.4 
GHz processor. This clearly demonstrates the high effi- 
ciency of FLMC simulations (with multiple occupancy of 
lattice sites and Kronecker (5-function interactions) over 
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both conventional lattice MC simulations (with SMAW 
and nearest-neighbor interactions) and FOMC simula- 
tions (with pair-potential calculations) . When compared 
with the corresponding lattice field theories based on the 
same Hamiltonian, FLMC simulations further provide a 
powerful means for unambiguously and quantitatively re- 
vealing the correlation/fluctuation effects in the system. 
Although here we use confined compressible homopoly- 



mers as a simple example, FLMC simulations can be 
readily performed for more complex systems such as poly- 
mer blends and block copolymers or even incompressible 
systems (where p(r) = po is required at all lattice sites 
but not SMAW). These will be the topics of our future 
publications. 

Financial support for this work was provided by NSF 
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